/*
 *  sparsematrixumf.cpp
 *  EPPI-0.0
 *
 *  Created by Juan Carlos García Orden on 07/07/2010.
 *
 */

#include "sparsematrixumf.h"

sparseMatrixUMF::sparseMatrixUMF ()
:	sparseMatrix("umf"),
	Symbolic(),
	Numeric(),
	b_out(),
	b_aux()
{
	setTypeName("umf");
}

// ***

sparseMatrixUMF::sparseMatrixUMF (int nrow, int ncol)
:	sparseMatrix("umf", nrow, ncol),
	Symbolic(),
	Numeric(),
	b_out(),
	b_aux()
{
	setTypeName("umf");
}

// ***

sparseMatrixUMF::~sparseMatrixUMF ()
{
	umfpack_di_free_symbolic(&Symbolic) ;
	umfpack_di_free_numeric(&Numeric);
}

// ***

void sparseMatrixUMF::initMatrices(void)
{
	b_out.resize(n);
//	b_outA.set_size(n);
	b_outAE.resize(n);
	b_outAEV.resize(n);
//	vaux.set_size(n);
	vauxE.resize(n);
}
	
// ***

void sparseMatrixUMF::factorize(int facMode)
{
	if(facMode==0)
	{
		(void) umfpack_di_symbolic(n,n,ja,ia,values,&Symbolic,NULL,NULL) ;
	}

	umfpack_di_free_numeric(&Numeric);
	
	(void) umfpack_di_numeric(ja,ia,values,Symbolic,&Numeric,NULL,NULL) ;	
}

// ***

void sparseMatrixUMF::subsSolve(vector<double>& b_in)
{
/*
	int i;
	cout << "ia is " << endl;
	for (i=0; i<nnz; ++i)	cout << ia[i] << "\t";
	cout << endl;
	
	cout << "ja is " << endl;
	for (i=0; i<n+1; ++i)	cout << ja[i] << "\t";
	cout << endl;
	
	cout << "values is " << endl;
	for (i=0; i<nnz; ++i)	cout << values[i] << "\t";
	cout << endl;
		
	cout << "b_in es " << endl;
	printSTLVector(b_in);
*/
	//if (lTwoNormSq(b_in)>=1.0e-6)
	//{
		(void) umfpack_di_solve(UMFPACK_A,ja,ia,values,&(b_out[0]),&(b_in[0]),Numeric,NULL,NULL) ;
	//}
	//else setZero(b_out);

	/*
	cout << "b_in es " << endl;
	printSTLVector(b_in);
	cout << "Norma es " << lTwoNormSq(b_in) << endl;

	cout << "b_out es " << endl;
	printSTLVector(b_out);
	*/
	
	b_in = b_out;	
}

// ***
//Arma
/*
void sparseMatrixUMF::subsSolve(rowvec & b_in)
{
	(void) umfpack_di_solve(UMFPACK_A,ja,ia,values,&(b_outA(0)),&(b_in(0)),Numeric,NULL,NULL) ;
	
	b_in = b_outA;
}
*/
// ***

void sparseMatrixUMF::subsSolve(RowVectorXd & b_in)
{
	(void) umfpack_di_solve(UMFPACK_A,ja,ia,values,&(b_outAE(0)),&(b_in(0)),Numeric,NULL,NULL) ;
	
	b_in = b_outAE;	
}

// ***

void sparseMatrixUMF::subsSolve(VectorXd & b_in)
{
	(void) umfpack_di_solve(UMFPACK_A,ja,ia,values,&(b_outAEV(0)),&(b_in(0)),Numeric,NULL,NULL) ;
	
	b_in = b_outAEV;	
}


// ***
/*
void sparseMatrixUMF::subsSolve(vector<vector<double> >& B_in)
{
	int i;
	int j;
	
	if (b_aux.size()==0)	b_aux.resize(n);
	
	for (j=0; j<B_in[0].size(); ++j)
	{
		for (i=0; i<n; ++i)	b_aux[i] = B_in[i][j];
		subsSolve(b_aux);
		for (i=0; i<n; ++i)	B_in[i][j] = b_aux[i];
	}
}
*/
// ***

void sparseMatrixUMF::subsSolveTrans(vector<vector<double> >& B_in)
{
	
	int i;
	#pragma omp parallel for private (i)
	for (i=0; i<B_in.size(); ++i)	subsSolve(B_in[i]);
	/*
	int i;
	for (i=0; i<B_in.size(); ++i)
	{
		(void) umfpack_di_solve(UMFPACK_A,ja,ia,values,&(b_out[0]),&(B_in[i][0]),Numeric,NULL,NULL) ;
		B_in[i] = b_out;
	}
	*/
}

// ***
/*
void sparseMatrixUMF::subsSolveTrans(double*)
{
}
*/
// ***
//Arma
/*
void sparseMatrixUMF::subsSolveTrans(mat & B_in)
{
	int i;
	for (i=0; i<B_in.n_rows; ++i)
	{
		vaux = B_in.row(i);
		(void) umfpack_di_solve(UMFPACK_A,ja,ia,values,&(b_outA(0)),&(vaux(0)),Numeric,NULL,NULL) ;
		B_in.row(i) = b_outA;
	}
}
*/
// ***

void sparseMatrixUMF::subsSolveTrans(MatrixXd & B_in)
{
	int i;
	for (i=0; i<B_in.rows(); ++i)
	{
		vauxE = B_in.row(i);
		(void) umfpack_di_solve(UMFPACK_A,ja,ia,values,&(b_outAE(0)),&(vauxE(0)),Numeric,NULL,NULL) ;
		B_in.row(i) = b_outAE;
	}	
}

// ***

void sparseMatrixUMF::initVectors(int nvec)
{}

// ***

void sparseMatrixUMF::copyMatrix( sparseMatrix *b)
{	
	int i;
	
	sparseMatrixUMF *aux = dynamic_cast<sparseMatrixUMF*> (b);

	aux->n = n;
	aux->nrhs = nrhs;
	aux->nnz = nnz;
	aux->ia = ia;
	aux->ja = ja;
//	for (i=0; i<nnz; ++i)	aux->ia[i] = ia[i];
//	for (i=0; i<n+1; ++i)	aux->ja[i] = ja[i];
		
	for (i=0; i<nnz; ++i)
	{
//		cout << "value[" << i << "] = " << values[i] << endl;
		aux->values[i] = values[i];
	}
	
	aux->Symbolic = Symbolic;
	aux->Numeric = Numeric;
	aux->b_out = b_out;
	aux->b_aux = b_aux;
	
	b = dynamic_cast<sparseMatrix*> (aux);
}

// ***
/*
void sparseMatrixUMF::freeN()
{
	umfpack_di_free_numeric(&Numeric);
}
*/